Image fusion scheme for differential phase contrast imaging

ABSTRACT

The latest progresses in breast imaging using differential phase contrast techniques pose the question of how to fuse multiple information sources, yielded by absorption, differential phase, and scattering signals, into a single, informative image for clinical diagnosis. It is proposed to use an image fusion scheme based on a multiple-resolution framework. The three signals are first transformed into multiple bands presenting information at different frequencies and then a two-step processing follows: (1) intra-band processing enhances the local signal-to-noise ratio using a novel noise estimation method and context modeling; and (2) inter-band processing weights each band by considering their characteristics and contributions, and suppressing the global noise level. The fused image, looking similar to a conventional mammogram but with significantly enhanced detail features, is reconstructed by inverse transform. The fused image is compatible with clinical settings and enables the radiologists to use their years of diagnosis experiences in mammography.

BACKGROUND OF THE INVENTION Field of the Invention

The present invention relates to a method and a system for differentialphase contrast imaging in the medical environment using an image fusionscheme.

Phase contrast and scattering-based X-ray imaging can access informationwhich are not detectable in conventional absorption-based methods, andtherefore are considered as promising tools to significantly improvemedical X-ray imaging. In the last few years, X-ray gratinginterferometry has been proved to be a suitable technique for both phasecontrast and scattering imaging. It can simultaneously measure theabsorption contrast (AC) image, the differential phase contrast (DPC)image and the dark-field (scattering) contrast (DFC) images (also calledsmall angle scattering contrast image), providing much richerinformation of the underlying sample than conventional X-rayradiography. This technique has been first successfully implemented onsynchrotron facilities and further on conventional X-ray tubes, pavingthe road for its broad usage in medical applications.

Among radiological applications, mammography is one of the techniquesthat may profit significantly from this technique because of theadvantages of phase contrast in imaging soft tissue at lower doses andthe additional diagnostic values that could be gained from thescattering signals. Recently, a research team domiziled at the PSIpresented the first investigation of native, non-fixed whole breastsamples using a grating interferometer based on a conventional X-raytube configuration (MammoDPC). The preliminary results confirmed thatthis novel technique provides clinically relevant information tocomplement and improve the diagnosis for mammography.

To evaluate this technology and further apply it in clinicalcircumstances, a rising question has been how to effectively fuse theAC, DPC, and DFC signals into one single, but more informative image.This image fusion problem is crucial for radiologists to adopt theMammoDPC technique in clinical applications because in practice it cansignificantly reduce the time for diagnosis and the complexity tointerpret the three different physical signals and enable theradiologists to use their multiple years of breast cancer diagnosticexperience. Further, the fused image can often provide intuitive,additional clinical information that are not apparent in the singleimages, and therefore results in more accurate diagnosis and consequentbetter treatment. It has been shown that RGB or IHS colour-coded imagefusion method is an efficient and natural way to represent the threesignals and especially suitable for human vision [1]. However so far, itdoes not comply with conventional gray-level images in mammography whichradiologists and medical doctors are used to work with, as well as thehigh-resolution (grey-level only) clinical monitors.

BRIEF SUMMARY OF THE INVENTION

It is therefore the objective of the present invention to provide amethod and a system for an image fusion scheme in differential phasecontrast imaging which perfectly fits into the existing diagnosticenvironment and routines in today's radiologists work.

This objective is achieved according to the present invention by amethod and a system for calculated image fusion in medical X-rayimaging, preferably mammography, wherein the absorption data, thedifferential phase contrast data and the (small-angle) scattering dataof the images of the underlying sample are merged by a fusion algorithminto one single, grey-level image with enhanced details.

This fusion strategy merges therefore the three signals into one singlegrey-level image which has a similar appearance as the conventionalmammograms but with significantly enhanced detail features yielded bythe DPC and DFC signals. Moreover, to ensure accurate diagnosis, thefusion procedure will neither sabotage the resolution nor increase thenoise level. Dramatic differences in noise properties of the threesignals have been observed in grating interferometer, which aretherefore handled carefully in the fusion algorithm.

In a preferred embodiment of the present invention, the image fusionalgorithm is based on a multiple-resolution (MR) framework, where theoriginal image is decomposed into several sub-images containing theinformation of the original image at different spatial frequencies,preferably Laplacian, wavelets or similar transformations which carryout the decomposition. The decomposition therefore emphasizes thecontribution of the DPC and the DFC signals to the fused image since themost interesting contributions of the DPC and DFC signals materializesrather in the high frequency domain than the contribution of AC signalwhich contributes rather in the low frequency domain.

Preferably, absorption data AC, the differential phase contrast data DPCand the (small-angle) scattering data DFC are obtained from X-rayinvestigations based on grating-based interferometry,analyzer-crystal-based imaging, coded aperture imaging or any imagingtechnique that record these three images. It has to be noted that theterms absorption data, absorption image, absorption signal rather meanthe same, namely the intensities of the X-ray imaging for the absorptioncaused alteration of the incident x-ray beam. Accordingly, this alsoapplies to the differential phase contrast and the (small-angle)scattering.

As mentioned above, the fusion algorithm will focus on generally spokenthree main steps:

a) step 1: the absorption, differential phase contrast and (small-angle)scattering images are transformed into a MR domain (for instance,wavelet domain) consisting of multiple levels (s denotes the levelindex), each level containing several sub-bands (o denotes the bandindex);b) step 2: the sub-band images are processed and merged in the MRdomain; andc) step 3: the merged image is reconstructed by the inverse MRtransform.

An advantageous feature for the execution of step 2 can be achieved whenstep 2 includes either an intra-band, or inter-band or both processing,wherein the fused sub-band F_(l)({circumflex over (t)}) can be generallyexpressed by

${{F_{l}\left( \hat{t} \right)} = {\sum\limits_{X \in {\lbrack{{AC},{DPC},{DFC}}\rbrack}}\;{{wn}_{X,l} \cdot {wb}_{X,l} \cdot {w_{X,l}\left( \hat{t} \right)} \cdot {D_{X,l}\left( \hat{t} \right)}}}},$where X represents one of the possible image types: the absorptioncontrast (AC), the differential phase contrast (DPC) or the dark-fieldcontrast (DFC); D_(X,l)({circumflex over (t)}) represents theunprocessed coefficients of band l; l denotes the pair (s,o), indicatinga certain band at level s and of band o. oε{HL,LH,HH} when wavelettransform is used for the decomposition. {circumflex over (t)}=(m,n)denotes the 2D coordinates in the sub-band image, and wherein theintra-band processing, which is represented by w_(X,l)({circumflex over(t)})·D_(X,l)({circumflex over (t)}), assigns a weighting factorw_(X,l)({circumflex over (t)}) to each coefficient within band l withthe purpose to increase the local signal-to-noise ratio (SNR); andwherein the inter-band processing gives a global weighting factorwn_(X,l)·wb_(X,l) to each band with the purpose of selecting usefulinformation according to the image characteristics and constraining theglobal noise level; specifically, wb_(X,l) denotes the band selectionweighting factor and wn_(X,l) denotes the noise constraint weightingfactor.

Preferably, the weighting factor w_(X,l)({circumflex over (t)}) in theintra-band processing is given by

${{w_{X,l}\left( \hat{t} \right)} = {\phi\left( \frac{\sigma_{Sw}^{({s,o})}\left( \hat{t} \right)}{\sigma_{Nw}^{({s,o})}\left( \hat{t} \right)} \right)}},$where φ(x) is a linear function which normalizes x into a certain rangewherein the ratio

$\frac{\sigma_{Sw}^{({s,o})}\left( \hat{t} \right)}{\sigma_{Nw}^{({s,o})}\left( \hat{t} \right)}$is defined as the local SNR for D_(X,l)({circumflex over (t)}) andσ_(Sw) ^((s,o))({circumflex over (t)}) and σ_(Nw) ^((s,o))({circumflexover (t)}) denote the local estimated signal strength and the standardnoise variance, respectively; for simplicity, the superscript (s,o) ofσ_(Sw) ^((s,o))({circumflex over (t)}) and σ_(Nw) ^((s,o))({circumflexover (t)}) is omitted; σ_(Sw)({circumflex over (t)}) andσ_(Nw)({circumflex over (t)}) are used in the following claims; thethree images types and their sub-bands are treated in the same way, sothe subscript X,l is also dropped.

In the fusion algorithm, the noise has to be watched attentively sincean incorrect handling of the noise could have the potential to corruptthe entire method. Therefore, as a further preferred embodiment of thepresent invention, the local standard noise variance σ_(Nw)({circumflexover (t)}) in the MR domain is estimated using the prior knowledge ofthe spatial pixel-wise noise variance σ_(Ns)({circumflex over (t)}) inthe original image:

${\left\lbrack {\sigma_{Nw}^{({s,o})}\left( \hat{t} \right)} \right\rbrack^{2} \approx {\frac{1}{2\pi}{\int_{- \pi}^{\pi}{{{{H^{({s,o})}\left( \hat{\omega} \right)}}^{2} \cdot {{A\left( \hat{\omega} \right)}}^{2}}e^{j\hat{\omega}\hat{t}}\ d\hat{\omega}}}}},$where H^((s,o))(ω) is the frequency response of the cascaded waveletfilters at scale s and orientation o, for example H^((s,o))(ω)=Π_(i-0)^(s-1)H(2^(i) ω)G(2^(s) ω), where G(ω) and H(ω) are the scaling andwavelet filters, respectively wherein the standard wavelet notation hasbeen used:A({circumflex over (ω)}) is the Discrete-time Fourier transform (DTFT)of σ_(Ns)({circumflex over (t)}), that isA({circumflex over (ω)})=∫σ_(Ns)({circumflex over (t)})e^(−j{circumflex over (ω)}{circumflex over (t)}) d{circumflex over (t)};σ_(Ns)({circumflex over (t)}) is measured and/or calculated in advance.

Preferably, the local signal strength σ_(Sw)({circumflex over (t)})represents the weight of the coefficient D({circumflex over (t)});

σ_(Sw)({circumflex over (t)}) is the local l_(p) norm of the N neighborcoefficients of D({circumflex over (t)}), which is

${{\sigma_{Sw}\left( \hat{t} \right)} = \left( {\sum\limits_{\hat{i} \subseteq N}\;\left\lbrack {D\left( \hat{i} \right)} \right\rbrack^{p}} \right)^{1/p}},$where D(î) may include D({circumflex over (t)}) itself or neighborcoefficients from its parent band (s−1,o) or son band (s+1,o); adedicated estimation of σ_(Sw)({circumflex over (t)}) is achieved byContext Modeling, where σ_(Sw)({circumflex over (t)}) is determined asfollowing:

-   -   i) a context Z({circumflex over (t)}) is assigned to each        coefficient D({circumflex over (t)}), which indicates the        significance of the coefficient and is defined as a weighted        average of the absolute values of the N neighbors of        D({circumflex over (t)}):Z({circumflex over (t)})=w^(T)u_(t),        where w and u_(t) are the vectorized N×1 weighting factor and        the vectorized absolute values of the N neighbor coefficients;    -   ii) the weighting factor w is statistically obtained by        least-square (LS) estimation through the whole sub-band:        w=(U ^(T) U)⁻¹ U ^(T) |D|,    -   where U is a M×N matrix with each row being u_(t) ^(T) for all        {circumflex over (t)} and D is the M×1 vector containing all the        coefficients D({circumflex over (t)}). M is the number of        coefficients in the sub-band;    -   after Z({circumflex over (t)}) is decided for each coefficient,        coefficients whose contexts Z({circumflex over (t)}) are close        in value are located and sorted. Let B_(t) denotes the set of        coefficients {D({circumflex over (t)})} whose context falls in a        similar value range. Those coefficients are considered to have        similar significances and their standard variance σ_(Sw) gives a        measurement of the significance,

${{\sigma_{Sw}\left( \hat{t} \right)} = \sqrt{\max\left( {{{\frac{1}{L}{\sum\limits_{B_{t}}\;{D\left( \hat{t} \right)}^{2}}} - {\sigma_{Nw}^{2}\left( \hat{t} \right)}},0} \right)}},$

-   -   where L is the number of the coefficients in B_(t).

Further, the band selection weighting factors is

${{{wb}_{{AC},l}(s)} = \frac{1}{1 + {\eta_{AC}\left( \frac{s}{S} \right)}^{2}}},{{{wb}_{{DPC},l}(s)} = \left\{ {\begin{matrix}{\frac{{\eta_{DPC}\left( \frac{s - s_{0}}{S} \right)}^{2}}{1 + \left( \frac{s - s_{0}}{S} \right)^{2}},} & {s > s_{0}} \\{0,} & {s \leq s_{0}}\end{matrix},{{{wb}_{{DFC},l}(s)} = \left\{ {\begin{matrix}{\frac{{\eta_{DFC}\left( \frac{s - s_{0}}{S} \right)}^{2}}{1 + \left( \frac{s - s_{0}}{S} \right)^{2}},} & {s > s_{0}} \\{0,} & {s \leq s_{0}}\end{matrix},} \right.}} \right.}$where S is the maximal decomposition level and s=0, 1, . . . , S is thelevel index; for multiple sub-bands within a certain level, the sameweighting factor is used; η_(AC), η_(DPC) and η_(DFC) are the image typerelated constants controlling the values of the weighting factors andtherefore weighting the contributions of each image type to the fusionimage; and s₀ represents a threshold frequency.

Furthermore, the global noise constraint weighting factors are inverselyproportional to the average noise level of the sub-band; for instance,

$\begin{matrix}{{wn}_{{AC},l} = \frac{1/{\overset{\_}{\sigma}}_{{AC},l}}{{1/{\overset{\_}{\sigma}}_{{AC},l}} + {1/{\overset{\_}{\sigma}}_{{DPC},l}} + {1/{\overset{\_}{\sigma}}_{{DFC},l}}}} \\{{wn}_{{DPC},l} = {\frac{{\overset{\_}{\sigma}}_{{AC},l}}{{\overset{\_}{\sigma}}_{{DPC},l}} \times {wn}_{{AC},l}}} \\{{wn}_{{DFC},l} = {\frac{{\overset{\_}{\sigma}}_{{AC},l}}{{\overset{\_}{\sigma}}_{{DFC},l}} \times {wn}_{{AC},l}}}\end{matrix},$where σ _(AC,l) is the average noise variance of AC image at band l,

${{\overset{\_}{\sigma}}_{{AC},l} = \sqrt{\frac{1}{N}{\sum\;{\sigma_{Nw}^{2}\left( \hat{t} \right)}}}};$andfor the DPC and DFC image, σ _(DPC,l) and σ _(DFC,l) are decided in thesame way.According to a further preferred embodiment of the present invention,the absorption data, the differential phase contrast data andsmall-angle scattering data were obtained from an arrangement forX-rays, in particular hard X-rays, for obtaining quantitative X-rayimages from a sample, comprising:

-   -   a. an X-ray source;    -   b. three gratings or at least two gratings;    -   c. a position-sensitive detector with spatially modulated        detection sensitivity having a number of individual pixels;    -   d. means for recording the images of the detector;    -   e. means for evaluating the intensities for each pixel in a        series of images, in order to identify the characteristics of        the object for each individual pixel as an absorption dominated        pixel and/or a differential phase contrast dominated pixel        and/or an X-ray scattering dominated pixel and for executing the        fusion algorithm;    -   f. wherein the series of images is collected by continuously or        stepwise rotating from 0 to π or 2π either the sample or the        arrangement and the source relative to the sample.

In this context of the afore-mentioned embodiment, the arrangement forX-rays can take the images either according to the so-called “near fieldregime” or according to the “Talbot-regime”, wherein preferably for thenear-field-regime operation, the distance between the gratings (G1, G2)is chosen freely within the regime, and for the Talbot-regime thedistance is chosen according to

$\mspace{20mu}{D_{n,{sph}} = {\frac{L \cdot D_{n}}{L - D_{n}} = \frac{{L \cdot n \cdot {p_{1}^{2}/2}}\eta^{2}\lambda}{L - {{n \cdot {p_{1}^{2}/2}}\eta^{2}\lambda}}}}$  where  n = 1, 3, 5  …  , and $\eta = \left\{ {\begin{matrix}1 & {{{if}\mspace{14mu}{the}\mspace{14mu}{phase}\mspace{14mu}{shift}\mspace{14mu}{of}\mspace{14mu} G_{1}\mspace{14mu}{is}\mspace{14mu}\left( {{2\; l} - 1} \right)\frac{\pi}{2}},} & {p_{2} = {\frac{L + D_{n,{sph}}}{L}p_{1}}} \\2 & {{{if}\mspace{14mu}{the}\mspace{14mu}{phase}\mspace{14mu}{shift}\mspace{14mu}{of}\mspace{14mu} G_{1}\mspace{14mu}{is}\mspace{14mu}\left( {{2\; l} - 1} \right)\pi},} & {p_{2} = {\frac{L + D_{n,{sph}}}{L}\frac{p_{1}}{2}}}\end{matrix},} \right.$where l=1, 2, 3, D_(n) is an odd fractional Talbot distance when theparallel X-ray beam is used, while D_(n,sph) is that when the fan orcone X-ray beam is used, L is the distance between the source (X-ray)and the first grating (G1).

Typically, in the grating structure of the arrangement for X-rays afirst granting is a line grating either as an absorption grating or aphase grating, wherein the phase grating is a low absorption grating butgenerating a considerable X-ray phase shift, the latter preferably of πor odd multiples thereof. Additionally, a second grating (G2) is a linegrating having a high X-ray absorption contrast with its period beingthe same as that of the self image of the first grating; the secondgrating being placed closely in front of the detector, or integrated init, with its lines parallel to those of the first grating.

The phase stepping with this arrangement for X-rays can be preferablyperformed by the shift of one grating with respect to the othersgratings in a direction perpendicular to the orientation of the gratinglines of the gratings.

For performing a tomographic reconstruction of the images, the fusionprotocol in the fusion algorithm is applied to multiple angularprojections. The presented fusion protocol is applied to tomographicallyreconstructed slices

BRIEF DESCRIPTION OF THE SEVERAL VIEWS OF THE DRAWING

Preferred embodiment of the present invention are hereinafter describedin more detail with reference to the attached drawings wherein:

FIG. 1 shows schematically the set-up of a typical grating-based x-rayinterferometer;

FIG. 2 shows a MammoDPC images for an ROI close to a carcinoma; (a) isthe AC image, (b) is the DFC image and (c) is the derivative of the ACimage, for fair comparison with the DPC image; (d) represents the DPCimage;

FIG. 3 illustrates the principle of the MR image fusion framework;

FIG. 4 shows quantitatively the band selection weighting factors used inthe image fusion algorithm;

FIG. 5 shows the comparison of the fusion image and the AC image of apatient with a large tumor mass and spiculations; (a) is the AC image ofthe whole breast; (b) is the corresponding fusion image; (c) are thedetails of the ROI shown in the AC image; and (d) is the same ROI of thefusion image in details;

FIG. 6 shows the comparison of the fusion image and the AC image of apatient with irregular parenchyma; (a) is the AC image of the wholebreast; (b) is the corresponding fusion image; (c) show the details ofthe ROI shown in the AC image; and (d) shows the same ROI of the fusionimage in details.

DESCRIPTION OF THE INVENTION

According to the present invention, a method and a system for calculatedimage fusion in medical imaging is disclosed. Presently, an image fusionscheme for differential phase contrast mammography based on themultiple-resolution (MR) method is hereinafter discussed in more detail.The framework using MR method for image fusion has been established byG. Piella [2]. It follows a “decomposition—processing andfusing—reconstruction” style. However for a particular application, likeMammoDPC here, principles for the ‘processing and fusing’ should beexplored and designed. The proposed fusion scheme handles the threesignals differently with particular care on adding more “signals” than“noises”. The AC, DPC, and DFC signals are first transformed intomultiple levels and bands, which represent image information atdifferent scales and frequencies, using wavelets decomposition. Then atwo-step processing follows: First, an intra-band processing enhancesthe local signal-to-noise ratio (SNR) and features using a novel noiseestimation method to model the noise and the context modeling to modelthe signal. In the second step, an inter-band processing weights eachband by considering their characteristics and contributions, as well assuppressing the global noise level. High frequency components of DPC andDFC images are preferred according to prior diagnostic knowledge whilethe information of the AC image is well kept. The resulting fusion imagehas similar appearance as conventional mammogram but significantlyenhanced detail features, and the noise level is maintained at the sametime.

Fundamentals of MammoDPC Method

The principle of grating interferometry for differential phase contrastimaging has been well established in the past few years. Brieflyspeaking, the system is usually composed of an X-ray tube, athree-grating interferometry, and a 2D detector as shown in FIG. 1. FIG.1 shows a two-gratings set-up (top) and three-gratings set-up (bottom)for X-ray imaging as for example disclosed in the European Patentapplications EP 2 400 891 A1 and EP 1 887 936 A1 which are herewithincorporated by reference.

Image Formation

During the data acquisition, one of the gratings is moved perpendicularto the grating lines within at least one gratings period (phase steppingapproach). For each pixel on the detector, a quasi-sine intensity curveis recorded. The AC, DPC, and DFC signals are obtained by an informationretrieving algorithm operating on phase stepping curves with and withoutthe object.

-   -   The phase stepping curve can be approximately expressed by

$\begin{matrix}{{{I_{Det}\left( {x,y} \right)} \approx {{a_{0}\left( {x,y} \right)} + {{a_{1}\left( {x,y} \right)}{\sin\left\lbrack {{2\pi\frac{n}{N_{ps}}} + {\Phi\left( {x,y} \right)}} \right\rbrack}}}},} & (1)\end{matrix}$where N_(ps) is the number of phase steps, n is the index of the step.

-   -   Then the AC, DPC, and DFC images are defined as:

$\begin{matrix}{{{AC} = {{{- \ln}\; T} = {- {\ln\left( \frac{a_{0,{obj}}}{a_{0,{bg}}} \right)}}}},{{DPC} = {\frac{1}{\pi}\left( {\Phi_{obj} - \Phi_{bg}} \right)}},{{DFC} = {{{- \ln}\; V} = {- {{\ln\left( {\frac{a_{1,{obj}}}{a_{0,{obj}}} \cdot \frac{a_{0,{bg}}}{a_{1,{bg}}}} \right)}.}}}}} & (2)\end{matrix}$

Signals from the scans with and without object (background scan) arelabeled by “obj” and “bg”, respectively. It follows that ranges of therecorded signals are ACε[0 ∞], DPCε[−1 1], DFCε [0 ∞]. T and v are thedefined as the transmission image and visibility degradation image whichwill be used in next section.

Noise Properties

Although the three signals are obtained simultaneously, their intrinsicnoise properties are different, which means that they cannot be treatedequally in the image fusion process. With a fine-tuned system, thequantum noise on the detector is considered to be the main noise source.In this case, the spatial pixel-wise noise variances could be calculatedand measured based on error propagation:

$\begin{matrix}{\begin{matrix}{{\sigma_{AC}^{2} = {\frac{f^{r}}{N_{ps}a_{o}^{r}}\left( {1 + {\frac{f^{s}}{f^{r}} \cdot \frac{1}{\overset{\_}{T}}}} \right){\overset{\_}{T}}^{2}}},} \\{{\sigma_{DPC}^{2} = {\frac{f^{r}}{2\pi^{2}v^{2}N_{ps}a_{o}^{r}}\left( {1 + {\frac{f^{s}}{f^{r}} \cdot \frac{1}{\overset{\_}{T}{\overset{\_}{V}}^{2}}}} \right)}},} \\{{\sigma_{DFC}^{2} = {{\frac{f^{r}}{v^{2}N_{ps}a_{o}^{r}}\left\lbrack {{v^{2}\left( {1 + {\frac{f^{s}}{f^{r}} \cdot \frac{1}{\overset{\_}{T}}}} \right)} + {2\left( {1 + {\frac{f^{s}}{f^{r}} \cdot \frac{1}{\overset{\_}{T}{\overset{\_}{V}}^{2}}}} \right)}} \right\rbrack}{\overset{\_}{V}}^{2}}},}\end{matrix},} & (3)\end{matrix}$where T and V are the mean intensity of the transmission image andvisibility degradation image for one pixel as defined in Eq. (2),respectively. a_(o) ^(r) is the mean intensity of the phase steppingcurve without object. v is the visibility at that pixel which can bemeasured with a background scan. f^(r) and f^(s) are noise relatedcoefficients with and without object, which could be measured Error!Reference source not found. v, f^(r) and f^(s) are known constants for atuned system.

The noise variances of DPC and DFC images are generally larger than theabsorption image. Considering the possible existences of other noisesources, this noise estimation is considered to be optimistic. However,it could provide useful prior knowledge for post-processing especiallybecause the spatial distribution of the noise is taken into account. Inorder to use this knowledge in the MR image fusion scheme, one needs totransform the spatial noise variances into the MR domain which will bediscussed later.

Additional Diagnostic Values

Image fusion schemes for MammoDPC should be driven by the additionaldiagnostic values from the DPC and DFC images. For the DPC image, themain gain is given by the high spatial frequency components of itbecause the differential phase signal is generally sensitive to theedges at tissues interfaces. The DFC image, on the other hand,contributes greatly in enhancing microcalcifications and fibrousconnective tissue visualization. An example is illustrated in FIG. 2.Both DPC and DFC images provide additional and complementary informationto the AC image. Besides, DPC image has the potential to distinguish atype of malignant conversion and its extension within normal breasttissue. The advantages of DFC image in imaging microcalcifications andfibrous tissue can play an important role in early breast cancerdetection because both microcalcifications and fibrous tissue arerelated to the development of malignant lesions. From the point of viewof image processing, those detail features (edges, microcalcificationsand fibrous structures) belong to the high frequency components of theimages.

Therefore, in the proposed image fusion scheme, it is aimed at takingpreferable high frequency information from DPC and DFC, whereas lowfrequency information is mostly taken from AC.

FIG. 2 shows a MammoDPC images for an ROI close to a carcinoma. Part (a)shows the AC image; part (b) is the DFC image; and part (c) is thederivative of the AC image, for fair comparison with the DPC image. Part(d) represents the DPC image. The fibrous tissues which are hardly seenin the AC image is highly enhanced, visualized and distinguished in theDFC image as indicated by the arrows and circles in (a) and (b). The DPCimage is sensitive to the edges and also provided additional detailfeatures which can not be seen in the AC image as demonstrated in thered rectangular in (c) and (d).

In the following the image fusion scheme for MammoDPC is explained inmore detail. The principle of the multiple resolution (MR) framework isillustrated in FIG. 3. Generally speaking, the framework includes threesteps:

Step 1: The images are first transformed into the MR domain, consistingof multiple levels. Each level might contain several sub-bands dependingon the decomposition methods, e.g. Laplacian decomposition or wavelettransform. Each sub-band contains the information of the original imageat a different spatial frequency.Step 2: The band images are then processed and merged within thetransformed domain according to the peculiarities of the application.Step 3: The fusion image is reconstructed by inverse MR transform.

The present scheme therefore mainly acts at the level of Step 2,including an intra- and an inter-band processing, while using a standardwavelet transform at Step 1 and 3.

Many decomposition methods have been studied and are used in variousapplications. Different types of wavelets, with different properties,could have been used, however this is not a subject for the discussionin this patent application. Instead, a general fusion scheme is proposedfor the particular application here, namely, differential phase contrastmammography.

Another advantage of the MR method worth mentioning is that bothdenoising and contrast enhancement, two common processing steps inradiography, can be applied easily within the same framework, providingthe possibility to further improve the fusion result according to therequirements set by the radiologists.

As mentioned earlier, FIG. 3 illustrates the principle of the MR imagefusion framework. Images are first decomposed into multiple bands(corresponding to image information at different resolution) and thosebands are processed and merged according to the applications. The fusionimage is finally reconstructed by inverse transform.

Intra-Band Processing

The intra-band processing is designed to increase local SNR (signalnoise resolution). A weighting factor is assigned to each waveletcoefficient:I _(AC,l)({circumflex over (t)})=w _(AC,SNR)({circumflex over (t)})×D_(AC,l)({circumflex over (t)})I _(DPC,l)({circumflex over (t)})=w _(DPC,SNR)({circumflex over (t)})×D_(DPC,l)({circumflex over (t)})I _(DFC,l)({circumflex over (t)})=w _(DFC,SNR)({circumflex over (t)})×D_(DFC,l)({circumflex over (t)})  (4)where D_(AC,l)({circumflex over (t)}) represents the originalcoefficients of band l in the wavelet domain for the AC image. l denotesthe pair (s,o), indicating a certain band at scale s and orientationoε{HL,LH,HH}. {circumflex over (t)}=(m,n) denotes the 2D coordinates ofa given coefficient in wavelet domain. The same notation rules apply toDPC and DFC images.

Here, the DPC image is taken as an example. The weighting factor isgiven by:

$\begin{matrix}{{{w_{{DPC},{SNR}}\left( \hat{t} \right)} = {\phi\left( \frac{\sigma_{Sw}\left( \hat{t} \right)}{\sigma_{Nw}\left( \hat{t} \right)} \right)}},} & (5)\end{matrix}$where φ(x) is a linear normalized function which normalizes the ratio

$\frac{\sigma_{Sw}\left( \hat{t} \right)}{\sigma_{Nw}\left( \hat{t} \right)}$into the range for a given band. Too large w_(DPC,SNR)({circumflex over(t)}) is observed to cause obvious artifacts in the fusion image. Tokeep the information from the AC image, it's not necessary to processD_(AC,l)({circumflex over (t)}). Therefore w_(AC,SNR)({circumflex over(t)}) could be fixed to 1.

The goal now is how to estimate the “signal term”σ_(Sw)({circumflex over(t)}) and the “noise term” σ_(Nw)({circumflex over (t)}) in the waveletdomain.

Noise Estimation

In order to use the spatial noise information (Eq. (3)) within the MRframework, the expressions of the noise variances in the wavelet domainneed to be deducted. Without loss of generality, an additive noise modelis assumed with known variance for each pixel, with the noises not beingspatially correlated and having local stationary property. The noise atposition {circumflex over (t)} can be modeled by a standard, whiteGaussian noise process ε({circumflex over (t)}):n({circumflex over (t)})=σ_(Ns)({circumflex over (t)})·ε({circumflexover (t)}),  (6)where σ_(Ns)({circumflex over (t)}) is the known standard noise variancein the spatial domain as expressed by Eq. (3). n({circumflex over (t)})is a real-value, zero mean random process. Note A({circumflex over (ω)})as the Discrete-time Fourier transform (DTFT) of σ_(s)({circumflex over(t)}),

$\begin{matrix}{{\sigma_{Ns}\left( \hat{t} \right)} = {\frac{1}{2\pi}{\int_{- \pi}^{\pi}{{A\left( \hat{\omega} \right)}e^{j\hat{\omega}\hat{t}}\ d{\hat{\omega}.}}}}} & (7)\end{matrix}$

Equation (6) now turns into

$\begin{matrix}{{{n\left( \hat{t} \right)} = {\frac{1}{2\pi}{\int_{- \pi}^{\pi}{{A\left( \hat{\omega} \right)}e^{j\hat{\omega}\hat{t}}\ d\;{F\left( \hat{\omega} \right)}}}}},} & (8)\end{matrix}$with F({circumflex over (ω)}) being the DTFT of ε({circumflex over(t)}).

The spatial noise variance could be therefore expressed by

$\begin{matrix}\begin{matrix}{\left\lbrack {\sigma_{Ns}\left( \hat{t} \right)} \right\rbrack^{2} = {E\left\lbrack {{n\left( \hat{t} \right)} \cdot {n\left( \hat{t} \right)}} \right\rbrack}} \\{= {\left( \frac{1}{2\pi} \right)^{2}{\int_{- \pi}^{\pi}{\int_{- \pi}^{\pi}{{E\left\lbrack \ {d\;{F\left( {\hat{\omega}}_{1} \right)}\ d\;{F\left( {\hat{\omega}}_{2} \right)}} \right\rbrack}{A\left( {\hat{\omega}}_{1} \right)}{A\left( {\hat{\omega}}_{2} \right)}e^{j{\hat{\omega}}_{1}\hat{t}}e^{j{\hat{\omega}}_{2}\hat{t}}}}}}} \\{= {\left( \frac{1}{2\pi} \right)^{2}{\int_{- \pi}^{\pi}{\int_{- \pi}^{\pi}{{E\left\lbrack \ {d\;{F\left( {\hat{\omega}}_{1} \right)}\ d\;{F\left( {\hat{\omega}}_{2} \right)}} \right\rbrack}{A^{*}\left( {\hat{\omega}}_{1} \right)}{A\left( {\hat{\omega}}_{2} \right)}e^{{- j}{\hat{\omega}}_{1}\hat{t}}e^{j{\hat{\omega}}_{2}\hat{t}}}}}}} \\{= {\frac{1}{2\pi}{\int_{- \pi}^{\pi}{{A^{*}\left( \hat{\omega} \right)}{A\left( \hat{\omega} \right)}e^{j\hat{\omega}\hat{t}}\ d\hat{\omega}}}}} \\{= {\frac{1}{2\pi}{\int_{- \pi}^{\pi}{{{A\left( \hat{\omega} \right)}}^{2}e^{j\hat{\omega}\hat{t}}\ d\hat{\omega}}}}}\end{matrix} & (9)\end{matrix}$

Here, the fact is used that n({circumflex over (t)}) is a real-valueprocess, A({circumflex over (ω)})=A*(−{circumflex over (ω)}) and thesign of {circumflex over (ω)}₁ is changed for integration. Also,E[dF({circumflex over (ω)}₁)dF({circumflex over (ω)}₂)]=2πδ({circumflexover (ω)}₁−{circumflex over (ω)}₂)d{circumflex over (ω)}₁d{circumflexover (ω)}₂ because F({circumflex over (ω)}) has orthogonal incrementsdue to the white noise process.

Following the same procedure, let H^((s,o))(ω) be the frequency responseof the cascaded wavelet filters at scale s and orientation o, forexample H^((s,o))(ω)=Π_(i-0) ^(s-1)H(2^(i) ω)G(2^(s) ω), where G(ω) andH(ω) are the scaling and wavelet filters, respectively. The waveletdomain noise variance at scale s and orientation o is approximatelygiven by:

$\begin{matrix}\begin{matrix}{\left\lbrack {\sigma_{Nw}^{({s,o})}\left( \hat{t} \right)} \right\rbrack^{2} = {E\left\lbrack {{n_{w}\left( \hat{t} \right)} \cdot {n_{w}\left( \hat{t} \right)}} \right\rbrack}} \\{\approx {\frac{1}{2\;\pi}{\int_{- \pi}^{\pi}{{{{H^{({s,o})}\left( \hat{\omega} \right)}}^{2} \cdot {{A\left( \hat{\omega} \right)}}^{2}}e^{j\;\hat{\omega}\;\hat{t}}d\hat{\omega}}}}}\end{matrix} & (10)\end{matrix}$

The connection between the spatial noise variance and the noise variancein the wavelet domain is setup by Eq. (7) and Eq. (10). For simplicity,for now on, the superscript (s,o) is omitted, σ_(Nw)({circumflex over(t)}) is used in the following sections.

Signal Estimation

To optimally fuse the DPC and DFC images into the AC image, theessential step is to extract all the informative features from them.Ideally, those features need to be decided by diagnostic practices.However, at this stage, a focus is laid on keeping important features bymodeling, which can be considered as “blind signal estimation”.

The signal is estimated using context modeling (CM), which is a powerfultool and has been effectively used in data compressing and imagedenoising. It takes into account the naturally, spatially varyingcharacteristics of an image. The fundamental idea is to group pixels ofsimilar nature but not necessarily spatially adjacent in the waveletdomain, and then gather statistical information from those pixels. Inour MammoDPC, the images have an average size larger than 5000×4000pixels, therefore the statistics can be considered to be robust.

For each coefficient D({circumflex over (t)}) in the wavelet domain(treating the three images in the same way, so the subscript is droppedin the following), a context Z({circumflex over (t)}) is assigned to it,which indicates the significance of the coefficient:Z({circumflex over (t)})=w ^(T) u _(t),  (11)where Z({circumflex over (t)}) is defined as a weighted average of theabsolute values of its p neighbors, and w and u_(t) are the vectorizedp×1 weighting factor and the absolute values of the neighborcoefficients.

The weighting factor w is obtained by least-square (LS) estimation inthe sub-band:w=(U ^(T) U)⁻¹ U ^(T) |D|,  (12)where U is a N×p matrix with each row being u_(t) ^(T) for all{circumflex over (t)} and D is the N×1 vector containing all thecoefficients D({circumflex over (t)}). N is the number of coefficientsin the sub-band. In our implementation, nine neighbors are used, eightof which are the neighborhood pixels in the same band, and the other oneis from the parent band.

After Z({circumflex over (t)}) is decided for each coefficient, thevariance of the random variable D({circumflex over (t)}) is estimatedfrom other coefficients whose context variables are close in value toZ({circumflex over (t)}). Coefficients whose contexts Z({circumflex over(t)}) are similar are located by applying a moving window on the sortedZ({circumflex over (t)}). Let B_(t) denote the set of points{D({circumflex over (t)})} whose context falls in the moving windows.Those coefficients are considered to have similar significances andtheir standard variance σ_(Sw) is calculated:

$\begin{matrix}{{{\sigma_{Sw}\left( \hat{t} \right)} = \sqrt{\max\left( {{{\frac{1}{L}{\sum\limits_{B_{i}}{D\left( \hat{t} \right)}^{2}}} - {\sigma_{Nw}^{2}\left( \hat{t} \right)}},0} \right)}},} & (13)\end{matrix}$where L is the number of the coefficients in the moving window andσ_(Nw) ²({circumflex over (t)}) is the pixel-wise noise variancediscussed above.Inter-Band Processing

In this procedure, instead of weighting each coefficient, a globalweighting factor is assigned to each band with the purpose to constrainthe global noise level and select useful information according to theimage characteristics.

-   -   The weighting scheme is given by        F _(l) =wn _(AC,l) ×wb _(AC,l) ×I _(AC,l) +wn _(DPC,l) ×wb        _(DPC,l) ×I _(DPC,l) +wn _(DFC,l) ×wb _(DFC,l) ×I        _(DFC,l),  (14)        where wb denotes the band selection weighting factor and wn        denotes the noise constraint weighting factor and will be        introduced in the following section. I_(AC,l), I_(DPC,l) and        I_(DFC,l) are the band images after the intra-band processing.        Band Selection

As mentioned above, the higher frequency components of DPC and DFCimages are preferred by the present fusion scheme according to theinvention. To make sure that the fusion image looks similar to the ACimage, the “zero” frequency (LL band) and their low frequency componentsare not used. This ensures that the histogram of the fusion image doesnot change dramatically from that of the AC image.

To weight the bands, the ideal weighting factors should be selected byevaluating the amount of useful information contained in the sub-band.This requires practical methods to extract clinical-relative features.We initially apply a weighting scheme which is a heuristicgeneralization of our Fourier-domain fusion method and extend it to fusethe DFC image. Low frequency components of DPC and DFC images are notused or suppressed while high frequency components of them arepreferred.

The band selection weighting factors are given by

$\begin{matrix}{{{{wb}_{A\; C}(s)} = \frac{1}{1 + {\eta_{A\; C}\left( \frac{s}{S} \right)}^{2}}},} & (15) \\{{{wb}_{DPC}(s)} = \left\{ \begin{matrix}{\frac{{\eta_{DPC}\left( \frac{s - s_{0}}{S} \right)}^{2}}{1 + \left( \frac{s - s_{0}}{S} \right)^{2}},} & {{s > s_{0}},} \\{0,} & {s \leq s_{0}}\end{matrix} \right.} & (16) \\{{{wb}_{DFC}(s)} = \left\{ \begin{matrix}{\frac{{\eta_{DFC}\left( \frac{s - s_{0}}{S} \right)}^{2}}{1 + \left( \frac{s - s_{0}}{S} \right)^{2}},} & {{s > s_{0}},} \\{0,} & {s \leq s_{0}}\end{matrix} \right.} & (17)\end{matrix}$where S is the maximal decomposition level, and S=10 in our case. s=0,1, . . . , S is the level index. For multiple sub-bands within a certainlevel, the same weighting factor is used.

η_(AC), η_(DPC) and η_(DFC) are the image type related constantscontrolling the values of the weighting factors and therefore theyweight the contributions of each image type to the fusion image. It isclear that the factors for AC image are decreasing with s increasing,while those of DPC and DFC are increasing, resulting in higher frequencycomponents to be preferred. s₀ represents the cutting frequency, whichmeans that levels below s₀ are not used in the fusion. This applies toDPC and DFC images. Smaller s₀ means more low frequency information fromDPC and DFC are involved in the fusion image, however it has been foundto affect the appearance of the fusion image significantly and to makethe fusion image un-interpretable by the radiologists.

η_(AC)=0 means all the information from the AC image is kept, whileη_(AC)>0 gives preferences to the low frequency components of AC image.The loss of the high frequency components could be compensated bymerging the high frequency components of the DPC image and helpssuppress the total noise Error! Reference source not found. Eq. (16) andEq. (17) are similar, however η_(DPC)>η_(DFC) would mean DPC image ispreferred to DFC image, and vice versa. In this paper, the contributionsof the DPC image are preferred because the edge enhancement effectsindeed improved the visual detection of many features at the level ofdiagnosis. Further, the DFC image contributes most to the noise in thethree signals, therefore n_(DFC) is chosen to be smaller than η_(DPC),and accordingly will influence the final, fusion image less.

It is noted here that the strategy for choosing the band weightingfactors could be very different according to different clinicalscenarios and is not limited by Eq. (15)-(17). In practice, multipleweighting factor tables could be pre-stored, selected and used togenerate the fusion image “on-the-fly”, which gives the possibilities totune the fusion image in diagnostic practices.

Global Noise Constraint

During the development of the algorithm it was observed that the noiselevel of the fusion image is mainly affected by the weighting on thebands and not by the pixel-wise weighting of the intra-band processing.To suppress the noise or at least keep its level similar to that of theAC image, a band-level weighting factor is assigned to each band.

For each band, the same noise estimation method described above isadopted. The average noise variance σ _(AC,l) of AC image for band l isgiven by

${\overset{\_}{\sigma}}_{{A\; C},l} = {\sqrt{\frac{1}{N}{\sum{\sigma_{Nw}^{2}\left( \hat{t} \right)}}}.}$The same rule applies for DPC and DFC images.

The noise-constraint weighting factors are then given by

$\begin{matrix}{{{wn}_{{A\; C},l} = \frac{1/{\overset{\_}{\sigma}}_{{A\; C},l}}{{1/{\overset{\_}{\sigma}}_{{A\; C},l}} + {1/{\overset{\_}{\sigma}}_{{DPC},l}} + {1/{\overset{\_}{\sigma}}_{{DFC},l}}}}{{wn}_{{DPC},l} = {\frac{{\overset{\_}{\sigma}}_{{A\; C},l}}{{\overset{\_}{\sigma}}_{{DPC},l}} \times {wn}_{{A\; C},l}}}{{{wn}_{{DFC},l} = {\frac{{\overset{\_}{\sigma}}_{{A\; C},l}}{{\overset{\_}{\sigma}}_{{DFC},l}} \times {wn}_{{A\; C},l}}},}} & (18)\end{matrix}$for AC, DPC, DFC images, respectively.

Intuitively, the weighting factor is inverse proportional to the averagenoise level of the sub-band. By doing this, the noise level in thefusion image can be suppressed.

The image fusion scheme according to the present invention was evaluatedby experimental MammoDPC data. The data was taken with a Talbot-Laugrating interferometer installed at Paul Scherrer Institute in Villigen,Switzerland. The interferometer was operated at the 5th Talbot distance,at a tube voltage of 40 kVp with mean energy of 28 keV and at a currentof 25 mA.

Image Fusion Parameters

The detector noise related coefficients were measured, where f^(r)=0.33and f^(s)=0.27. The spatial noise variances of the three signals werecalculated according to Eq. (3). Decimated wavelets decomposition withDaubechies filters with order 4 was used. The size of three images was7000(H)×5000(W). The images were decomposed into 10 levels, and eachlevel contained three bands (LH,HL,HH) except the low frequency levelLL.

The band selection weighting factors were decided according to Eq.(15)-(17), where η_(AC)=0.5, η_(DPC)=2, η_(DFC)=1, and s₀=2 were used.These weighting factors are shown in FIG. 4. Band No. 0 represents theLL band (the “zero” frequency component).

A fusion example using the proposed scheme is given in FIG. 5. Thebreast of this patient contains a large tumor mass and manyspiculations. As it can be seen, the fusion image shows an overallbetter performance than the AC image in terms of sharpness and detailvisibility. The fusion scheme clearly enhances the spiculations whichare more obviously seen in the ROIs shown by FIGS. 5(c) and 5(d).Diagnostically, this enhancement is very relevant since spiculations areusually appearing in relation to malignant formations sometimesdifficult to detect in normal, absorption based mammogram. Some airbubbles artefacts are more obvious in the fusion image, however this arerelated to the intrinsically ex-vivo character of the experiment andwill disappear once in-vivo experiment will be performed. Therefore,removing those artifacts is not an issue and will not be discussedfurther.

FIG. 5 shows the comparison of the fusion image and the AC image of apatient with a large tumor mass and spiculations. (a) The AC image ofthe whole breast; (b) The corresponding fusion image; (c) The details ofthe ROI shown in the AC image; (d) The same ROI of the fusion image indetails. The images are shown with the same display window. Note alsothat both images have been acquired with the same dose.

Evaluation of an image fusion method is a tough topic and there are nostandard metric and/or method available. The evaluation strongly dependson the applications themselves. Here, two ROIs (a uniform area withinthe sample and one in the background), estimated the noise propertiesand compared the fusion image and the AC image to evaluate the noiseproperty of the proposed fusion scheme. The two ROIs are indicated inFIG. 5(b). The mean and standard variance of the ROI within the sampleare 0.913±0.008 in the AC image and 0.913±0.010 in the fusion image,respectively. For the ROI in the background, the values are 0.038±0.111in the AC image, and 0.039±0.110 in the fusion. This rough qualitativeanalysis shows that the noise level in the fusion image was wellpreserved.

Another example is shown in FIG. 6. This patient shows an irregularparenchyma because of a previous bilateral breast reduction mammoplasty.The diagnostic challenge was to distinguish invasive tumor from scartissue. The microscopic structure of the scars was expected to generatea larger scattering signal than normal breast tissue. These enhancedtissue features contributed by the DFC images (as indicated by thearrows in FIG. 6(b)) are effectively merged into the fusion image, whichprovide complementary and useful information for diagnosis. FIG. 6 showsthe comparison of the fusion image and the AC image of a patient withirregular parenchyma. (a) The AC image of the whole breast; (b) Thecorresponding fusion image. (c) The details of the ROI shown in the ACimage; (d) The same ROI of the fusion image in details. The enhanceddetails in the fusion image are pointed out by arrows. The images areshown with the same display window.

CONCLUSION

The effective image fusion scheme for grating-based differential phasecontrast mammography (MammoDPC) with the MR framework according to thepresent invention has been demonstrated. Higher frequencies of the DPCand DFC images are efficiently merged into the AC image, while noise isconstrained. Future work will include the image quality evaluation withclinical practices and fusion parameter optimization for differentapplications.

REFERENCES

-   [1] M. Stampanoni, Z. Wang, T. Thuering et al., The first analysis    and clinical evaluation of native breast tissue using differential    phase-contrast mammography, Invest. Radiol. 46 (2011) 801.-   [2] G. Piella, A general framework for multiresolution image fusion:    from pixels to regions, Information fusion 4 (2003) 259.-   [3] C. Clavijo, Z. Wang and M. Stampanoni, Wavelet-based noise model    driven denoising algorithm for differential phase contrast    mammography, Opt. Express 21 (2013) 10572.-   [4] S. G. Chang, B. Yu and M. Vetterli, Spatially adaptive wavelet    thresholding with context modeling for image denoising, IEEE Trans.    Image Process. 9 (2002) 1522.

The invention claimed is:
 1. A method for calculated image fusion inmedical X-ray imaging, which comprises the steps of: merging absorptiondata, differential phase contrast data and small-angle scattering dataof images of an underlying sample by a fusion algorithm into one single,grey-level image with enhanced details; basing the calculated imagefusion on a multiple-resolution framework, where an original image isdecomposed into several sub-images containing information of theoriginal image at different spatial frequencies; obtaining theabsorption data, the differential phase contrast data and thesmall-angle scattering data of the images from X-ray investigationsbased on grating-based interferometry, analyzer-crystal-based imaging,or coded aperture imaging, and using an arrangement for X-rays forobtaining quantitative X-ray images from a sample, the arrangement forX-rays containing: an X-ray source; at least two gratings; and aposition-sensitive detector with spatially modulated detectionsensitivity having a number of individual pixels; transforming theabsorption data, the differential phase contrast data and thesmall-angle scattering data of the images into a multiple-resolution(MR) domain consisting of multiple levels (s denotes a level index),each level containing several sub-bands (o denotes a band index);processing and merging sub-band images in the MR domain by performingintra-band processing, inter-band processing or both the intra-bandprocessing and the inter-band processing, wherein a fused sub-bandF_(l)({circumflex over (t)}) can be generally expressed by${{F_{l}\left( \hat{t} \right)} = {\sum\limits_{X \in {\lbrack{{A\; C},{DPC},{DFC}}\rbrack}}{{wn}_{X,l} \cdot {wb}_{X,l} \cdot {w_{X,l}\left( \hat{t} \right)} \cdot {D_{X,l}\left( \hat{t} \right)}}}},$where X represents one of possible image types: an absorption contrast(AC), a differential phase contrast (DPC) or a dark-field contrast(DFC); D_(X,l)({circumflex over (t)}) represents an unprocessedcoefficients of band l; l denotes a pair (s,o), indicating a certainband at level S and of band o, oε{HL,LH,HH} when a wavelet transform isused for a decomposition; {circumflex over (t)}=(m,n) denotes 2Dcoordinates in a sub-band image: the intra-band processing, which isrepresented by w_(X,l)({circumflex over (t)})·D_(X,l)({circumflex over(t)}), assigns a weighting factor w_(X,l)({circumflex over (t)}) to eachcoefficient within the band l to increase a local signal-to-noise ratio;the inter-band processing gives a global weighting factorwn_(X,l)·wb_(X,l) to each band for selecting useful informationaccording to image characteristics and constraining a global noiselevel; and wb_(X,l) denotes a band selection weighting factor andwn_(X,l) denotes a noise constraint weighting factor; and reconstructinga merged image by an inverse multiple-resolution transform.
 2. Themethod according to claim 1, wherein the weighting factorw_(X,l)({circumflex over (t)}) in the intra-band processing is possiblygiven by${{w_{X,l}\left( \hat{t} \right)} = {\phi\left( \frac{\sigma_{Sw}^{({s,o})}\left( \hat{t} \right)}{\sigma_{Nw}^{({s,o})}\left( \hat{t} \right)} \right)}},$where φ(x) is a linear function which normalizes x into a certain rangewherein a ratio$\frac{\sigma_{Sw}^{({s,o})}\left( \hat{t} \right)}{\sigma_{Nw}^{({s,o})}\left( \hat{t} \right)}$is defined as the local SNR for D_(X,l)({circumflex over (t)}) andσ_(Sw) ^((s,o))({circumflex over (t)}) and σ_(Nw) ^((s,o))({circumflexover (t)}) denote a local estimated signal strength and a local standardnoise variance, respectively; for simplicity, superscript (s,o) ofσ_(Sw) ^((s,o))({circumflex over (t)}) and σ_(Nw) ^((s,o))({circumflexover (t)}) is omitted; σ_(Sw)({circumflex over (t)}) andσ_(Nw)({circumflex over (t)}) are used in the following claims; thethree images types and their sub-bands are treated in a same way, sosubscript X,l is also dropped.
 3. The method according to claim 2,wherein the local standard noise variance σ_(Nw)({circumflex over (t)})in the MR domain is estimated using prior knowledge of a spatialpixel-wise noise variance σ_(Ns)({circumflex over (t)}) in the originalimage:${\left\lbrack {\sigma_{Nw}^{({s,o})}\left( \hat{t} \right)} \right\rbrack^{2} \approx {\frac{1}{2\;\pi}{\int_{- \pi}^{\pi}{{{{H^{({s,o})}\left( \hat{\omega} \right)}}^{2} \cdot {{A\left( \hat{\omega} \right)}}^{2}}e^{j\;\hat{\omega}\;\hat{t}}d\hat{\omega}}}}},$where H^((s,o))(ω) is a frequency response of cascaded wavelet filtersat scale s and orientation o, for example H^((s,o))(ω)=Π_(i-0)^(s-1)H(2^(i) ω)G(2^(s) ω), where G(ω) and H(ω) are scaling and waveletfilters, respectively wherein the standard wavelet notation has beenused: A({circumflex over (ω)}) is the Discrete-time Fourier transform(DTFT) of σ_(Ns)({circumflex over (t)}), that isA({circumflex over (ω)})=∫σ_(Ns)({circumflex over (t)})e^(−j{circumflex over (ω)}{circumflex over (t)}) d{circumflex over (t)};and σ_(Ns)({circumflex over (t)}) is measured and/or calculated inadvance.
 4. The method according to claim 2, wherein: the localestimated signal strength σ_(Sw) ^((s,o))({circumflex over (t)})represents a weight of the coefficient D({circumflex over (t)});σ_(Sw)({circumflex over (t)}) is the local l_(p) norm of the N neighborcoefficients of D({circumflex over (t)}), which is${{\sigma_{Sw}\left( \hat{t} \right)} = \left( {\sum\limits_{\hat{i} \subseteq N}\left\lbrack {D\left( \hat{i} \right)} \right\rbrack^{p}} \right)^{1/p}},$where D(î) may include D({circumflex over (t)}) itself or neighborcoefficients from its parent band (s−1,o) or son band (s+1,o); adedicated estimation of σ_(Sw)({circumflex over (t)}) is achieved byContext Modeling, where σ_(Sw)({circumflex over (t)}) is determined asfollowing: (i) a context Z({circumflex over (t)}) is assigned to eachsaid coefficient D({circumflex over (t)}), which indicates asignificance of the coefficient and is defined as a weighted average ofabsolute values of the N neighbor coefficients of D({circumflex over(t)}):Z({circumflex over (t)})=w ^(T) u _(t), where w and u_(t) are avectorized N×1 weighting factor and vectorized absolute values of the Nneighbor coefficients; (ii) The weighting factor w is statisticallyobtained by least-square estimation through a whole sub-band:w=(U ^(T) U)⁻¹ U ^(T) |D|, where U is a M×N matrix with each row beingu_(t) ^(T) for all {circumflex over (t)} and D is the M×1 vectorcontaining all the coefficients D({circumflex over (t)}), M is a numberof the coefficients in the sub-band; (iii) After Z({circumflex over(t)}) is decided for each coefficient, coefficients whose contextsZ({circumflex over (t)}) are close in value are located and sorted; LetB_(t) denotes a set of coefficients {D({circumflex over (t)})} whosecontext falls in a similar value range, those coefficients areconsidered to have similar significances and their standard varianceσ_(Sw) gives a measurement of the significance,${{\sigma_{Sw}\left( \hat{t} \right)} = \sqrt{\max\left( {{{\frac{1}{L}{\sum\limits_{B_{t}}{D\left( \hat{t} \right)}^{2}}} - {\sigma_{Nw}^{2}\left( \hat{t} \right)}},0} \right)}},$where L is the number of the coefficients in B_(t).
 5. The methodaccording to claim 1, wherein the band selection weighting factors is${{{wb}_{{A\; C},l}(s)} = \frac{1}{1 + {\eta_{A\; C}\left( \frac{s}{S} \right)}^{2}}},{{{wb}_{{DPC},l}(s)} = \left\{ {{\begin{matrix}{\frac{{\eta_{DPC}\left( \frac{s - s_{0}}{S} \right)}^{2}}{1 + \left( \frac{s - s_{0}}{S} \right)^{2}},} & {{s > s_{0}},} \\{0,} & {s \leq s_{0}}\end{matrix}{{wb}_{{DFC},l}(s)}} = \left\{ \begin{matrix}{\frac{{\eta_{DFC}\left( \frac{s - s_{0}}{S} \right)}^{2}}{1 + \left( \frac{s - s_{0}}{S} \right)^{2}},} & {{s > s_{0}},} \\{0,} & {s \leq s_{0}}\end{matrix} \right.} \right.}$ where S is a maximal decomposition leveland s=0, 1, . . . , S is a level index; for multiple sub-bands within acertain level, a same weighting factor is used; η_(AC), η_(DPC) andη_(DFC) are an image type related constants controlling values of theweighting factors and therefore weighting contributions of each imagetype to the fusion image; and s₀ represents a threshold frequency. 6.The method according to claim 1, wherein global noise constraintweighting factors are inversely proportional to an average noise levelof a sub-band; for instance,${wn}_{{A\; C},l} = \frac{1/{\overset{\_}{\sigma}}_{{A\; C},l}}{{1/{\overset{\_}{\sigma}}_{{A\; C},l}} + {1/{\overset{\_}{\sigma}}_{{DPC},l}} + {1/{\overset{\_}{\sigma}}_{{DFC},l}}}$${wn}_{{DPC},l} = {\frac{{\overset{\_}{\sigma}}_{{A\; C},l}}{{\overset{\_}{\sigma}}_{{DPC},l}} \times {wn}_{{A\; C},l}}$${{wn}_{{DFC},l} = {\frac{{\overset{\_}{\sigma}}_{{A\; C},l}}{{\overset{\_}{\sigma}}_{{DFC},l}} \times {wn}_{{A\; C},l}}},$where σ _(AC,l) is an average noise variance of AC image at band l,${{\overset{\_}{\sigma}}_{{A\; C},l} = \sqrt{\frac{1}{N}{\sum{\sigma_{Nw}^{2}\left( \hat{t} \right)}}}};$and for a DPC and DFC image, σ _(DPC,l) and σ _(DFC,l) are decided in asame way.
 7. The method according to claim 1, wherein the arrangementfor X-rays further comprising: a processing device to record the imagesof the position-sensitive detector; the processing device evaluatingintensities for each pixel in a series of images, in order to identifycharacteristics of an object for each individual pixel as one of anabsorption dominated pixel, a differential phase contrast dominatedpixel, or an X-ray scattering dominated pixel; and wherein a series ofimages is collected by continuously or stepwise rotating from 0 to π or2π either the sample or the arrangement and the x-ray source relative tothe sample.
 8. The method according to claim 7, wherein the at least twogratings include a first grating being a line grating selected from thegroup consisting of an absorption grating and a phase grating, whereinthe phase grating is a low absorption grating but generating aconsiderable X-ray phase shift, the phase grating being of π or oddmultiples thereof.
 9. The method according to claim 8, wherein: the atleast two gratings include a second grating being a line grating havinga high X-ray absorption contrast with its period being a same as that ofa self image of the first grating; and the second grating being placedclosely in front of the position-sensitive detector, or integrated init, with its lines parallel to those of the first grating.
 10. Themethod according to claim 9, which further comprises performing phasestepping by a shift of one of the gratings with respect to other ones ofthe gratings in a direction perpendicular to an orientation of gratinglines of the first and the second gratings.
 11. The method according toclaim 7, which further comprises applying a fusion protocol to multipleangular projections acquired with a purpose of performing a tomographicreconstruction.
 12. The method according to claim 7, wherein a presentedfusion protocol is applied to tomographically reconstructed slices. 13.The method according to claim 1, wherein the images are taken accordingto a so-called “near field regime” or according to a “Talbot-regime”.14. The method according to claim 13, wherein for a near-field-regimeoperation, a distance between the gratings is chosen freely within thenear field regime, and for the Talbot-regime is chosen according to$\mspace{20mu}{D_{n,{sph}} = {\frac{L \cdot D_{n}}{L - D_{n}} = \frac{{L \cdot n \cdot {p_{1}^{2}/2}}\;\eta^{2}\lambda}{L - {{n \cdot {p_{1}^{2}/2}}\;\eta^{2}\lambda}}}}$  where  n = 1, 3, 5  …  , and $\eta = \left\{ {\begin{matrix}1 & {{{if}\mspace{14mu}{the}\mspace{14mu}{phase}\mspace{14mu}{shift}\mspace{14mu}{of}\mspace{11mu} G_{1}\mspace{14mu}{is}\mspace{14mu}\left( {{2l} - 1} \right)\frac{\pi}{2}},} & {p_{2} = {\frac{L + D_{n,{sph}}}{L}p_{1}}} \\2 & {{{if}\mspace{14mu}{the}\mspace{14mu}{phase}\mspace{14mu}{shift}\mspace{14mu}{of}\mspace{14mu} G_{1}\mspace{14mu}{is}\mspace{14mu}\left( {{2l} - 1} \right)\pi},} & {p_{2} = {\frac{L + D_{n,{sph}}}{L}\frac{p_{1}}{2}}}\end{matrix},} \right.$ where l=1, 2, 3 . . . , D_(n) is an oddfractional Talbot distance when a parallel X-ray beam is used, whileD^(n,sph) is that when a fan or cone X-ray beam is used, L is a distancebetween the X-ray source and the first grating.
 15. The method accordingto claim 1, wherein the decomposing into the several sub-images isperformed by one of Laplacian, wavelets or similar transformations tocarry out the decomposition.